First-principles study of ferroelectricity and isotope effects in H-bonded KDP crystals 
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By means of extensive first-principles calculations we studied the ferroelectric phase transition and 
the associated isotope effect in KH2PO4 (KDP). Our calculations revealed that the spontaneous 
polarization of the ferroelectric phase is due to electronic charge redistributions and ionic displace- 
ments which are consequence of proton ordering, and not vice versa. The experimentally observed 
double-peaked proton distribution in the paraelectric phase cannot be explained by a dynamics of 
only protons. This requires, instead, collective displacements within clusters that include also the 
heavier ions. These tunneling clusters can explain the recent evidence of tunneling obtained from 
Compton scattering measurements. The sole effect of mass change upon deuteration is not sufficient 
to explain the huge isotope effect. Instead, we find that structural modifications deeply connected 
with the chemistry of the H-bonds produce a feedback effect on tunneling that strongly enhances 
the phenomenon. The resulting influence of the geometric changes on the isotope effect agrees with 
experimental data from neutron scattering. Calculations under pressure allowed us to analyze the 
issue of universality in the disappearance of ferroelectricity upon compression. Compressing DKDP 
so that the distance between the two peaks in the deuteron distribution is the same as for protons in 
KDP, corresponds to a modification of the underlying double-well potential, which becomes 23 meV 
shallower. This energy difference is what is required to modify the O-O distance in such a way as to 
have the same distribution for protons and deuterons. At the high pressures required experimentally, 
the above feedback mechanism is crucial to explain the magnitude of the geometrical effect. 



I. INTRODUCTION. 

Potassium dihydrogen phosphate (KH2PO4, or KDP) 
crystals are a key component in quantum electronics. 
They are widely used in controlling and modulating the 
frequency of laser radiation in optoelectronic devices, 
amongst other uses such as TV screens, electrooptic de- 
flector prisms, i nterdigital electrodes, light deflectors, 
and adjustable light filters. Besides the obvious techno- 
logical interest, KDP is also interesting from a fundamen- 
tal point of view. KDP is a prototype ferroelectric (FE) 
crystal belonging to the family of hydrogen-bonded fer- 
roelectrics, which was extensively studied in the past. 1 ' 2 
Their PO4 molecular units are linked by hydrogen bonds, 
and ferroelectricity appears to be connected to the behav- 
ior of the protons in these H-bonds. Normal (H2O) ice is 
the most prominent member of this family, 3,4 which also 
includes other compounds like PbHPCU, 5 and squaric 
acid (C4H2O4). 6 What makes KDP particularly interest- 
ing is the possibility of growing quite large, high-quality 
single crystals from solution, thus making it very suit- 
able for experimental studies. Indeed, a large wealth of 
experimental data has been accumulated during the sec- 
ond half of the past century. 1,2 

Phosphates in KDP are linked through approximately 
planar H-bonds forming a three-dimensional network. In 
the paraelectric (PE) phase at high temperature, the H- 
atoms occupy with equal probability two symmetrical po- 
sitions along the H-bond separated a distance S, which 
characterizes the so-called disordered phase. Below the 



critical temperature T c w 122 K, the protons localize 
into one of the symmetric sites, thus leading to the or- 
dered FE phase. Here, the spontaneous polarization P s 
appears perpendicular to the proton ordering plane, and 
the PO4 tetrahedra distort. The proton configuration in 
this phase is depicted in Fig. 1; each PO4 unit has two 
covalently bonded and two H-bonded hydrogen atoms, 
following Slater's ice rules 8 . The oxygen atoms that bind 
covalently to the hydrogen are called acceptor (02 in Fig. 
1), and those H-bonded are called donors (01 in Fig. 1). 

A striking feature common to all H-bonded ferro- 
electrics is undoubtedly the huge isotope effect observed 
upon deuteration. In fact, the deuterated compound 
(DKDP) exhibits a T c about two times larger than KDP. 
This giant effect was first explained by the quantum 
tunneling model proposed in the early sixties 9 . Within 
the assumption of interacting, single-proton double wells, 
this model proposes that individual protons tunnel be- 
tween the two wells. Protons are more delocalized than 
deuterons, thus favoring the onset of the disordered PE 
phase at a lower T c . Improvements of the above model 
include coupling between the proton and the K-PO4 dy- 
namics. 10-13 These models have been validated a poste- 
riori on the basis of their predictions, although there is 
no direct experimental evidence of tunneling. Only very 
recent neutron Compton scattering experiments seem to 
indicate the presence of tunneling 14 . However, the con- 
nection between tunneling and isotope effect remains un- 
clear, in spite of recent careful experiments. 15 

On the contrary, a series of experiments carried out 
since the late eighties 16-20 provided increasing experi- 



mental evidence that the geometrical modification of the 
hydrogen bonds and the lattice parameters upon deuter- 
ation (Ubbelohde effect 21 ) is intimately connected to the 
mechanism of the phase transition. The distance d be- 
tween the two collective equilibrium positions of the pro- 
tons (see Fig. 1) was shown to be remarkably correlated 
with T c . 19 Actually, it seems that proton and host cage 
are connected in a non-trivial way, and are not sepa- 
rable. 22 These findings stimulated new theoretical work 
where virtually the same phenomenology could be ex- 
plained without invoking tunneling. 23-25 However, these 
theories were developed at a rather phenomenological 
level. Only very recently, the first ab initio calculations, 
based on Density Functional Theory (DFT), were con- 
ducted in these systems. 26-30 These approaches have the 
advantage of allowing for a confident and parameter-free 
analysis of the microscopic changes affecting the different 
phases in this system. 

In this work we investigate, using DFT electronic 
structure calculations within the generalized gradient ap- 
proximation to exchange and correlation, the relationship 
between proton ordering, internal geometry, polarization, 
tunneling and isotope effects in KDP (details of the meth- 
ods used are exposed in Section II). To this end, consid- 
eration of the following questions naturally arises: (1) 
What is the microscopic mechanism which gives rise to 
the FE instability?, (2) How do local instabilities lead to 
the double-site distribution in the PE phase?, (3) What is 
the quantum origin of the geometrical effect?, (4) What 
is the main cause of the giant isotope effect: tunneling or 
the geometrical modification of the H-bonds?, (5) How 
does pressure affect the energetics and the structural pa- 
rameters in the system? 

With the aim of shedding light on the above formu- 
lated questions, and on the general problematic in KDP, 
we conducted different computational experiments and 
made a revision of previously obtained results. First, we 
carried out electronic structure calculations in the tetrag- 
onal unpolarized phase (PE), forcing protons to be in the 
middle of the O-H-0 bonds. Calculations in the polar- 
ized phase (FE), with the H ordered off-center, were per- 
formed in both, the tetragonal fixed cell and in the com- 
pletely relaxed cell, which is orthorhombic. We studied 
the structures and the charge reorganization leading to 
the FE instability. These results are presented in Section 
III. In Section IV, we analyze global instabilities in KDP 
to understand the relation between proton ordering and 
polarization. To address the tunneling issue, we stud- 
ied local instabilities by determining the dependence of 
the system energetics upon the proton position in the H- 
bonds under various conditions: allowing or not K and P 
ions relaxations, and considering also individual proton 
and small cluster displacements. Besides, we show in this 
section a calculation of the momentum distribution of the 
proton along the H-bond in different phases and compare 
the results with recent experimental data. Section V is 
devoted to a thorough study of quantum fluctuations, 
and the controversial problem of the isotope effects. Wc 



show how a self-consistent quantum modeling, based on 
our first principles calculations, is able to explain the 
striking mass dependence of the geometrical effect. In 
Section VI, we present calculations of the energetics and 
the structural parameters as a function of pressure. We 
show that the results of related experiments under pres- 
sure are explained by the non- linear relationship between 
deuteration and geometric effects, derived in the previ- 
ous section. Finally, in Section VII we discuss the above 
issues, and elaborate our conclusions. 

II. AB INITIO METHODS 

We have performed ab initio calculations of KDP, 
within the framework of DFT, 31:32 using two different 
pseudopotential codes, one based on localized basis sets 
(LB), and another using plane waves (PW). 

O 




FIG. 1. Schematic view of the internal structure of KDP 
along the tetragonal axis. The fractional coordinates of P and 
K atoms along the c axis, are indicated in brackets. Covalent 
and H-bonded hydrogens are connected to corresponding oxy- 
gens by full and broken lines, respectively. 

The LB calculations were carried out using the SIESTA 
program. 33,34 This is a fully self-consistent DFT method 
that employs a linear combination of pseudoatomic or- 
bitals (LCAO) of the Sankey-Niklewsky type as basis 
functions 35 . These basis functions are strictly confined 
in real space, what is achieved by imposing in the pseu- 
doatomic problem (i.e. the atomic problem where the 
Coulomb potential has been replaced by the same pseu- 
dopotential that will be used in the solid state), the 
boundary condition that the orbitals vanish at a finite 
cutoff radius, rather than at infinity as for the free 
atom. Therefore, these solutions are slightly different 
from the free atom case, and have somewhat larger as- 
sociated energies because of the confining potential. The 
relevant parameter for this approximation is precisely 
the orbital confinement energy E c , given by the energy 
difference between the eigenvalues of the confined and 
the free orbitals. In our calculations, we used a value 
of E c = 50mcV. By decreasing this value further, we 
checked that we obtain total energies and geometries with 
sufficient accuracy. For the representation of the valence 



TABLE I. Comparison of the ab-initio (LB and PW) internal structure parameters with DKDP experimental data 41 for 
the different cases considered in the text. The notation is the same used in the experimental works referred. 7 is the relative 
z-displacement of the K and P atoms from the equidistant situation (see definition in section IV. A) . Distances in A and angles 
in degrees. 
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electrons in the LB method we used double-zeta bases 
with polarization functions (DZP). This means two sets 
of orbitals for the angular momenta occupied in the 
isolated atom, and one set of orbitals for the first 
non-occupied angular momentum (polarization orbitals). 
Again, this size of the basis set turns out to be accurate 
enough for our purposes. 

The exchange-correlation energy terms were computed 
using the Perdew-Burke-Ernzerhof (PBE) form of the 
generalized gradient approximation 36 . This type of func- 
tionals has already been used to describe hydrogen- 
bonded systems, with quite good accuracy. 37 We also 
tried the BLYP functional 38 , which gives very good 
results for molecular systems. However, the results 
in the solid state were of quality inferior than PBE. 
We used non-local, norm-conserving Troullier-Martins 
pseudopotentials 39 to eliminate the core electrons from 
the description. We also included nonlinear core cor- 
rections (NLCC) for a proper description of the K ion, 
due to an important overlap of the core charge with the 
valence charge density in this atom. We checked this 
approximation by comparing with a K + pseudopoten- 
tial that includes the 3s and 3p shells explicitly in the 
valence, as semicore states (9 explicit electrons per K 
atom). The semicore results are very closely reproduced 
by the core-corrected calculations. For the real-space grid 
used to compute numerically the Coulomb and exchange- 
correlation integrals 33 ' 34 , we used an equivalent energy 
cutoff of 125 Ry. 

The pseudopotcntial PW calculations were carried out 
with the PWSCF code 40 , with the same exchange- 
correlation functional and pseudopotentials, except for 
the core-corrected K, where we used a semicore pseu- 
dopotential for K + (9 electrons per K atom). The plane 
wave expansion was cut off at a maximum PW kinetic 
energy of 150 Ry. Such a high cutoff was necessary to 



obtain convergence in the internal degrees of freedom, 
particularly the hydrogen-bonded units. The PWSCF 
code also allows for the computation, within the linear 
response regime, of vibrational and dielectric properties 
such as phonon frequencies, Born effective charges and 
dielectric constant. Phonon eigenvectors were used to 
calculate the total energy curves under pressure, by con- 
straining the optimization to motions preserving the pat- 
tern of the ferroelectric normal mode, which is related to 
the parameter S. 

The PE phase of KDP has a body-centered tetrago- 
nal (bet) structure with 2 formula units per lattice site 
(16 atoms). For the LB calculations that describe homo- 
geneous distortions, we used the conventional bet cell (4 
formula units), but doubled along the tetragonal c axis. 
This supercell comprises 8 formula units (64 atoms). A 
larger supercell is required to describe local distortions. 
To this end, we used the equivalent conventional fct cell 
(containing 8 formula units, and axes rotated through 45 
degrees with respect to the conventional bet cell), also 
doubled along the c-axis (128 atoms). The LB calcu- 
lations were conducted using a r-point sampling of the 
Brillouin zone (BZ). This choice of sampling proved suf- 
ficient provided the large supercells used in the calcula- 
tions. 

Most of the calculations have been done using the LB 
approach which, within the approximations described 
above, turned out into quite a fast computational pro- 
cedure, compared to PW calculations. As a test, we 
checked the LB approach against the PW results. The 
PW calculations were carried out on the 16-atom bet 
unit cell, with a BZ sampling consisting of 8 centered 
Monkhorst-Pack k-points. This number of points was 
checked for convergence, and proved sufficient. The re- 
sults for the geometrical parameters are reported in ta- 
ble I. It can be seen that the LB values are of quality 



comparable to the PW results. The differences can be 
attributed mainly to the approximation made with the 
confinement of the orbitals in the LB calculations. 



III. CHARACTERIZATION OF THE 
STRUCTURES AND CHARGE FLOW 
MECHANISMS 

We first optimized the structure with paraelectric 
phase symmetry. To this end, we constrained the H- 
atoms to remain centered in the O-H-0 bonds, and fixed 
the lattice parameters to the experimental values of the 
deuterated compound (DKDP) at T c +5 K (a = b = 7.459 
A and c = 6.957 A) in the conventional bet cell. 41 The 
choice of DKDP instead of KDP for the comparison with 
experiment is based on the fact that nuclear quantum ef- 
fects, which are neglected in the first-principles calcula- 
tions, are less important. Optimization of all the atomic 
positions leads to what we call the unpolarized tetragonal 
(UT) structure. This can be interpreted as an average of 
the true paraelectric phase. In fact, according to experi- 
mental data, in this latter the H-atoms are observed with 
equal probability in two symmetric off-centered positions 
along the H-bonds. In Table I we compare the relevant 
structural parameters resulting from both types of cal- 
culations and also experimental data. The agreement 
between the two theoretical approaches is quite good - 
thus validating the later use of the LB approach -, and 
their comparison with experiment is very satisfactory, ex- 
cept for the 0-0 distance doo which, specially in the UT 
case, turns out to be too short. 26 This delicate issue will 
be discussed below. 

Maintaining the lattice parameters and constraining 
the K and P atoms to their centered positions in the UT 
structure, we next allowed for H off-center relaxation to- 
wards the ordered configuration sketched in Fig. 1. The 
0-0 distance is also optimized. In this way we obtained 
a H off-center shift 5/2 = 0.154 A and an 0-0 distance 
of 2.472 A. We will show below in more detail that the 
H off-centering produces an electronic charge redistribu- 
tion from the neighborhood of the 02 atoms towards that 
of the 01 atoms. As a consequence, unbalanced forces 
are generated that favo pairing of the K and P atoms 
along the z axis, on the charge-excess side (01) of the 
PO 4 units. The former observation indicates that, con- 
straining the K and P to their centered positions, does 
not prevent the H atoms from abandoning the center of 
the H-bonds. The centered position for the H atoms is 
always unstable, as we will show in the next Section. 

The next step was to relax also the positions of the 
K and P atoms, thus leading to the polarized tetrag- 
onal (PT) structure, whose geometrical parameters are 
listed in Table I. This PT structure is not yet the ground 
state, because the ferroelectric distortion is coupled to a 
shear strain mode. It is this acoustic mode that becomes 
soft before the FE mode, thus piloting a structural tran- 
sition to an orthorhombic phase, which is very similar 



to the PT. This coupling can be observed in the a xyi 
off-diagonal components of the calculated stress tensor. 
According to this observation, we relaxed again all the 
internal degrees of freedom, but now fixing the simula- 
tion cell to the experimental orthorhombic structure of 
DKDP at T c — 10 K. This corresponds to lattice param- 
eters a = 10.598 A, b = 10.496 A and c = 6.961 A in 
the conventional fct cell. 41 The calculated geometrical 
parameters, which are close to those of the PT structure, 
are shown in the last three columns in Table I compared 
to experimental data. 

In the experimental orthorhombic structure, however, 
the stress tensor is diagonal but not isotropic. This indi- 
cates that, if the lattice parameters were also to be opti- 
mized, the b/a ratio would be different from the experi- 
mental value. In addition, the isotropic part of the stress 
(the pressure) is non-zero, thus indicating a small dif- 
ference in equilibrium volume between calculations and 
experiment. In general the agreement is quite reasonable, 
again with the exception of the 0-0 distance, which is 
0.1 A too small in the UT structure. This is a very im- 
portant issue, because the potential for the deuterons (or 
protons) in the H-bond is extremely sensitive to the 0-0 
distance 42 . In the present DFT-PBE calculations for the 
UT structure this distance is 2.42 A, i.e. 0.1 A shorter 
than the experimental value 41 . This difference, which 
can even change the shape of the potential felt by the 
deuteron in the H-bond, cannot be fully attributed to 
the optimization for centered deuterons. In fact, it per- 
sists when we optimize the structure in the orthorhombic 
FE phase, although slightly reduced (2.49 A vs. 2.53 A). 
One possible reason are quantum nuclear effects. Our 
calculations are for clamped nuclei, corresponding to in- 
finite deuteron mass. If quantum dynamics of deuterons 
was to be included, it would slightly increment this dis- 
crepancy because nuclear derealization favors shorter H- 
bonds. This can also be seen from the fact that the ex- 
perimental 0-0 distance for protons is shorter than for 
deuterons (Ubbelohde effect). Therefore, the inclusion of 
quantum effects would imply even shorter 0-0 distances. 
It is neither a problem of the pscudopotcntial approach, 
which has been tested against all-electron calculations. 

We conclude, then, that the main origin of the under- 
estimation of the 0-0 bond length is in the approximate 
character of the exchange-correlation functional. In fact, 
calculations for related gas- phase systems like H 3 0^ in- 
dicate a similar 0.06 A underestimation when compar- 
ing GGA values to correlated quantum chemical calcula- 
tions 43 . Moreover, present test calculations for the water 
dimer also indicate and underestimation of the doo dis- 
tance by 0.06 A with respect to experimental values 44 . 

Therefore, the differences in the H-bond geometry app- 



TABLE II. Changes q(PT) - q(UT) in the Mulliken orbital and bond overlap populations in going from the UT to the PT 
configuration, in units of e/1000. 
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FIG. 2. Differential charge density contours Ap(r) in the planes containing the following atoms: (a) P-Ol- • H, (b) P-02-H. 
Labels 01 and 02 denote the positions of the respective nuclei, positioned at (0,0). Labels 02-P and Ol-P indicate the position 
of the center of the corresponding bonds. The same convention is used for the 02-H and 01- ■ -H bonds. Positive (negative) 
contours are in solid (dashed) lines. The thickest lines represent an absolute value of 2.96 x 10 -3 eA -3 . The thinner lines are 
obtained by successively halving this value, down to 3.70 x 10~ 4 eA -3 . 



ear to be mainly due to the approximate character of the 
exchange-correlation term. Unfortunately, at present, a 
sufficiently well validated and efficient scheme to go be- 
yond GGA is lacking. Therefore, in order to avoid prob- 
lems derived from this feature in the study of the stability 
of local cluster distortions in next Section, we decided to 
fix the 0-0 distances in the host to the experimental 
values. 

With the purpose of analyzing the charge redistribu- 
tions produced by the ordered proton off-centering, we 
computed the changes in the Mulliken orbital and bond- 
overlap populations in going from the UT to the PT con- 
figuration, as shown on Table II. Mulliken populations 
depend strongly on the choice of the basis set. Differ- 
ences, however, are much less sensitive. 

An increase of the charge localized around 01 can be 
clearly observed; the main contribution (» 70%) is pro- 
vided by a decrease in the 02 charge. 26 The trends ob- 
served in Table II are confirmed by the charge density 
difference Ap(r) = ppt{y) — put{y). In Fig. 2a and 2b 
we plot cuts of the above quantity in the planes deter- 
mined by the atoms P-OT ■ H and P-02-H , respectively. 
A combined analysis of both, Table II and Fig. 2a, in- 



dicates a significant enhancement of the population of 
the 01 atom, accompanied by a smaller increment in the 
Ol-P orbitals. This happens at the expenses of the pop- 
ulation of the 01- • H and 02-P overlap orbitals, and the 
population of the 02 atom. Therefore, as two H-atoms 
move away from 01 and other two approach 02, the 01- 
H bond weakens and the 02-H bond strengthens. The 
charge localizes mostly around 01 and, to a lesser ex- 
tent, in the P-Ol orbitals. This is consistent with the 
increase of d(01 — 02) and the decrease of d(P — Ol) 
reported in Table I. The contrary occurs in the vicinity 
of the 02 atom, as indicated by the orbital and bond- 
overlap populations in Table II and the contours in Fig. 
2b. The overall effect is a flow of electronic charge from 
the 02 side of the tetrahedron towards the 01 side, and 
a concomitant modification of its internal geometry. The 
charge redistribution is rather local, and gives rise to a 
polarization composed by electronic and ionic contribu- 
tions . This polarization, whose origin can therefore be 
traced back to the off-centering of the H-atoms in the 
perpendicular plane, is intimately linked to ferroelectric- 
ity. In fact, the combined motion of all the atoms and 
the concomitant electronic redistribution corresponds to 



an unstable phonon in the UT structure . When this 
phonon mode freezes into one of the two stable minima, 
we obtain the PT structure which has a polarization, and 
is thus ferroelectric. A schematic view of this combined 
effect is presented in Fig. 3. 
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FIG. 3. Schematic view, perpendicular to the c-axis, of the 
atomic motions (solid arrows) and electronic charge redistri- 
butions (dotted curved arrows) , happening upon off-centering 
of the H-atoms. The percentages of the total charge redis- 
tributed are also shown for the charge-transfers occurring be- 
tween different orbitals and atoms. 

A further confirmation of these ideas comes from the 
saturated polarization which was obtained from linear 
response calculations of Born effective charges. 46 ' 40 Con- 
sidering the eigenvector of the FE mode e a (i) , where a 
indicates the Cartesian coordinate and i the ion, we cal- 
culated the dynamical effective charge in the z direction: 

z;(fe) = J2Y, z ^ 7= = L6 e - W 

i a * 

The effective charge components in the x and y directions 
vanish. Z* (FE) is multiplied by the FE-mode amplitude 
corresponding to the stable minimum, giving rise to a sat- 
uration polarization of P s = 3.25/iC/cm 2 , slightly lower 
than experimental values. Analyzing the individual con- 
tributions, we observe that a substantial part arises, in 
fact, from the P ions. However, this contribution is only a 
40% of the total polarization. The other 60% arises from 
the H atoms through a non-diagonal xz (or yz) compo- 
nent of the effective charge tensor. This component is 
Z* Z (H) = 0.6 e while Z* Z (P) = 3 e, but the displace- 
ment of the H atoms in the x (or y) axis is more than 
five times larger than that of the P ions along the z axis. 
This, and the fact that there are twice as many H than 
P in the unit cell, explains why H contributes as much 
as P to the polarization. The interesting observation is 
that the H atoms move in the plane perpendicular to 
the z axis. Therefore, their off-diagonal contribution can 
only be due to electronic polarization effects. Although 



the effective charges of the K and O atoms are not null, 
their small displacements lead to negligible contributions 
to the spontaneous polarization. 

IV. GLOBAL AND LOCAL FERROELECTRIC 
INSTABILITIES 

A. Correlation Between Proton Ordering and 
Polarization 

In the previous section we showed that there is an in- 
stability of the system towards the PT structure as hy- 
drogens collectively move away from the centers of the O- 
H-0 bond. From our simulations, we observe that when 
the protons are constrained to remain centered in the H- 
bonds, the K and P atoms are stable in their centered 
positions. However, centering the heavy atoms does not 
imply the centering of the H-atoms. Protons, in fact, are 
never stable at their centered positions. This provides 
a strong evidence that the origin of ferroelectricity is in 
the off-center ordering of the protons, and that proton 
off-centering and ferroelectricity are very correlated phe- 
nomena. 

To identify the driving mechanism of the ferroelectric 
instability, we analyzed the relationship between proton 
ordering and polarization. To this purpose, we inves- 
tigated the ab initio potential energy surface (PES) as 
a function of the proton off-centering parameter 6 = 
doo — 2dou, and the K-P relative displacement along 
the c-axis, which we quantify in terms of the parameter 
7 = dpp — 2oIkp, with oIkp the smallest K-P distance. 
It is worth mentioning here that 7 is a measure of po- 
larization, since a test calculation provided us a linear 
relationship between these quantities. 

We fully relax the oxygen positions for each chosen 
(<5, 7) pair, and plot the energy contours of the bidi- 
mensional PES in the inset to Fig. 4. The charac- 
teristics of this PES are as follows: it exhibits a sad- 
dle point at 6 = 7 = 0, and two equivalent minima at 
(5, 7) ~ ±(0.3, 0.15) A On the one hand, from the energy 
contours it can be seen that at 5 — (centered protons) 
there is no instability for any value of 7, i.e. the crystal 
is stable against polarization (7 7^ 0) unless the protons 
are ordered off-center. On the other hand, even for van- 
ishing 7 (polarization) the energy minimum corresponds 
to a finite S, i.e. protons are always collectively unsta- 
ble at the H-bond center. This is further visualized in 
Fig. 4, where we plot the energy profiles as a function of 
6 for different fixed values of 7. For 7 = 0, the energy 
profile exhibits a double-well in the S coordinate with a 
barrier of ~ 6 meV per molecular unit. For increasing 
values of 7 the minima are always at 6 ^ 0, up to a value 
of 7 w 0.02 A, where one of the two minima completely 
disappears. Therefore, we conclude from the above con- 
siderations that the source of the ferroelectric instability 
is the H off-centering, and not viceversa. 
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FIG. 4. Energy profiles as a function of 5, for values of 
7 = (solid line), 0.02 (short-dashed), 0.05 (dot-dashed) and 
0.1 (long-dashed) A. The inset shows equispaced energy con- 
tours (step = 13.6 meV/ KH2PO4 unit). The minima at 
(S, 7) ~ ±(0.3, 0.15)A lie ~ 50 meV below the saddle point at 
(0,0) 



B. Local Instabilities, Quantization in Small 
Clusters, and the Nature of the Paraelectric Phase 

We now address the microscopic origin of the observed 
proton double-occupancy in the PE phase 41 , which is an 
indication of the order-disorder character of the transi- 
tion. This phenomenon can be ascribed either to static 
or thermally activated dynamic disorder, or to tunneling 
between the two sites. Any of these possibilities requires 
the search for instabilities with respect to correlated but 
localized H motions in the PE phase, including also the 
possibility of heavy-ions relaxation. In fact, correlated 
motions of a large number of protons become increas- 
ingly unlikely in a tunneling scenario, because this im- 
plies higher barriers and heavier effective masses, thus 
reducing the tunneling probability. To analyze localized 
distortions we consider increasingly larger clusters em- 
bedded in a host paraelectric matrix. For the reasons 
exposed in the previous Section, the host is modeled by 
protons centered between oxygens, and the experimental 
structural parameters (including the 0-0 distances) of 
KDP at Tf DP +5 K (127 K). 41 In order to assess the 
effect of the volume increase observed upon deuteration, 
we also analyze the analogous case of D in DKDP by 
expanding the host structural parameters to the corre- 
sponding experimental values at T® KDP +5 K (234 K). 

41 

We analyze results for different clusters comprising N 
hydrogens (deuteriums): (a) N=l H(D) atom, (b) N=4 
H(D) atoms which connect a PO4 group to the host, 



(c) N=7 H(D) atoms localized around two PO4 groups, 
and (d) N=10 H(D) atoms localized around three PO4 
groups. For all these clusters we consider correlated mo- 
tions with the pattern shown in Fig. 1, which are the 
most favorable for exhibiting FE instabilities, as it was 
illustrated in the previous section. This correlated pat- 
tern, is represented by a single collective coordinate x 
whose value coincides with the H(D) off-center displace- 
ment 5/2. Two cases are considered: (i) first, we allow 
for the motion of H atoms alone, maintaining all other 
atoms fixed, (ii) second, we also allow for the relaxation 
of the heavy ions K and P, which follow the ferroelectric 
mode pattern 47,46 , as expected. Subsequent quantization 
of the cluster motion in the corresponding effective po- 
tential allows for the determination of the importance of 
tunneling in the disordered phase. Rigorously, the size 
dependence should be studied for larger clusters than 
those mentioned here. However, it will be shown be- 
low that short- range quantum fluctuations in the PE 
phase are sufficiently revealing, especially far away from 
the critical point. 

In Fig. 5 we show, for the clusters considered, the 
total energy variation as a function of x. For the case of 
H motions alone, we do not observe any instability for 
N=l and N=4, both in KDP and DKDP. A small barrier 
of ~ 6 meV appears in DKDP for the N=7 move, as 
shown in Fig. 5(b) (open squares). This barrier grows 
up to 25 meV for the N=10 cluster in DKDP (open 
circles) . 




x [A] x [A] 

FIG. 5. Energy profiles for correlated local distortions in 
(a) KDP and (b) DKDP. Reported are clusters of: 4 H(D) 
(diamonds), 7 H(D) (squares), and 10 H(D) (circles). Empty 
symbols and dashed lines indicate that only the H(D) atoms 
move. Motions that involve also heavy atoms (P and K) are 
represented by filled symbols and solid lines. Lines are guide 
to the eye only. 

However, quantum mechanical calculations of the cluster 



levels, which will be described bellow, yield ground states 
(GS) quantized above the barriers and, consequently, the 
absence of tunneling in those cases (see Fig. 5). In KDP, 
even the largest cluster considered is very stable, as in- 
dicated by the open circles in Fig. 5(a). The result for 
KDP suggests to rule out this type of motions in the 
paraelectric phase, because they are incompatible with 
double site occupancy. 

In a second step we considered also the motion of the 
heavy atoms in the above correlated local motions. The 
situation changed drastically, as shown by the solid lines 
and full symbols in Fig. 5(a) and (b). In fact, clusters in- 
volving two or more PO4 units - cases (c) and (d) above 
- exhibit instabilities in both KDP and DKDP, with a 
significant barrier in DKDP for case (d), of the order of 
w 150 meV. We note here that the instability appears 
in clusters which are sufficiently large, thus providing a 
measure of the FE correlation length. Moreover, the in- 
stabilities are much stronger (and the correlation length 
accordingly shorter) in the expanded DKDP lattice, than 
in KDP. 

We treat these clusters quantum-mechanically, by solv- 
ing the Schrodingcr equation for the collective coordinate 
x. This is done for each cluster in the corresponding 
effective potentials of Fig. 5. The effective mass for 
the local collective motion of the cluster is calculated 
as n = mid?, where i runs over the displacing atoms 
and rrii are their corresponding atomic masses, at is the 
i-atom displacement at the minimum from their posi- 
tions in the PE phase, relative to the H(D) displacement. 
The effective masses per H(D) calculated for these cor- 
related motions in different clusters are about /in ~ 2.3 
([Id ~ 3.0) proton masses (m p ) in KDP (DKDP), re- 
spectively The calculation of the GS energy in the heavy 
clusters, leads now to quantized levels below the barri- 
ers, as shown by dotted lines in Fig. 5(b), for all clusters 
in DKDP. This is a clear sign for tunneling arising from 
correlated D motions involving also heavy ions. These 
collective motions can be understood as a local distortion 
reminiscent of the global FE mode. 46 In KDP, however, 
even the largest cluster considered (N=10), has the GS 
level quantized above the barrier. The onset of tunneling 
at a critical cluster size, provides a rough indication of 
the correlation volume: it comprises more than 10 hydro- 
gens in KDP, but no more than 4 deuteriums in DKDP. 
49 We clearly observe, then, that the dynamics of the 
order-disorder transition would involve fairly large H(D)- 
clusters together with heavy-atom (P and K) displace- 
ments. Thus, the observed proton double-occupancy is 
explained in our calculation by the tunneling of large 
and heavy clusters. The last conclusion is confirmed by 
the double-site distribution determined experimentally 
for the P atoms. 50 ' 51 

A possible scenario for the FE phase transition would 
then be the following: the PE phase is made of clusters 
of different size, some of them (the large ones) preserve 
the FE structure, i.e. do not exhibit double site H(D) 
occupancy because the barrier is too high. Other clus- 



ters (smaller) have lower barriers and smaller effective 
masses, and thus can tunnel, giving rise to double oc- 
cupancy. The effect of temperature is that of modifying 
the preferential cluster size, which grows as a measure 
of the correlation length on approaching the transition. 
When the average cluster size reaches a value in which 
neither tunneling nor thermal hopping are anymore al- 
lowed, then the phase transition takes place. Of course, 
this is a mean field vision, but we believe that the picture 
is quite plausible. 

C. Momentum Distributions 

Since Blinc's model proposal, 9 it has been subject of 
controversy whether the protons are actually tunneling 
between the two equivalent sites along the bridges, or 
they are localized in one of these sites and jump to the 
other through phonon assisted tunneling in the paraelec- 
tric phase. Reiter et al. 14 have recently attempted to 
elucidate this question by performing neutron Compton 
scattering experiments. Due to the much shorter time 
scale of this experiment, compared to typical times for 
phonon assisted jumps in the paraelectric phase, it is 
claimed that it is possible to distinguish between a pro- 
ton coherently distributed between the two equivalent 
sites, and one which is alternatively occupying one site 
or the other. In this experiment, the momentum dis- 
tribution along the bridge, n(p), has been obtained by 
inverting the measured scattering function under plau- 
sible conditions. 14 Very significant changes in n(p) are 
observed when going through the transition, which were 
not to be expected if the proton was localized only in 
one of the equivalent sites, in both phases. As shown by 
the solid lines in Fig. 6, n(p) is considerably narrower 
in the high temperature phase, indicating an increase in 
the spread of the region where the proton is coherently 
distributed (the wave packet). More conclusively, the 
high temperature distribution shows a zero and a subse- 
quent oscillation which correspond precisely to a double 
peaked spatial wavefunction, i.e. the proton coherently 
distributed over both sites along the bond. In contrast, 
well below T c , n(p) shows a single and broader maximum 
at p = 0, thus indicating single-site occupancy. 

Our calculations for the coherent motion of hydro- 
gen clusters with fixed heavy ions, in a host of a mean 
paraelectric phase, indicate that only very large clusters 
would exhibit double well potentials with energy barri- 
ers high enough to allow for collective tunneling. On the 
other hand, considerably smaller clusters are able to tun- 
nel if also the heavy ions are allowed to move coherently 
with H in KDP, or with D in DKDP. Therefore, Compton 
scattering results can be explained if the observed coher- 
ent double peaked distribution of a proton along a bridge 
is interpreted as part of a coherent motion together with 
heavy ions in a cluster. 

Since the largest cluster we treated in KDP is not able 



to tunnel, a direct comparation of our momentum distri- 
bution calculations in the PE phase with experiment is 
not feasible. Instead, we can make a prediction of what 
would the n(p) distribution be like, in the PE phase of 
DKDP. For this purpose, we considered the correspond- 
ing double-well potentials as functions of the position 
of D along the bridge, for the 4-D and 7-D clusters in 
DKDP. We calculated the wavefunctions with the cluster 
effective masses [i — 10.4 m p and 21.4 m p , respectively. 
The resulting momentum distributions of DKDP in the 
PE phase are shown in Fig. 6 (right panel), together with 
the the experimental curve of KDP for illustrative pur- 
poses. The second oscillation arises from the quantum 
coherence of the real space distribution. As the effec- 
tive mass of the cluster increases, the second oscillation 
has larger amplitudes, while the main oscillation remains 
unchanged. 

In the ferroelectric phase of KDP the proton distri- 
bution is single peaked, corresponding to a single-well 
anharmonic potential for each proton. The momentum 
distribution calculated for a single proton wave function 
in such potential is shown by the dashed line in the left 
panel of Fig. 6. In order to understand the difference 
with the experimental data, we performed another cal- 
culation where the surrounding ions are allowed to re- 
lax. This leads to a shallower potential, but the effective 
mass also increases. The result, shown as a dotted line, 
deviates even more from the data. The deviation from 
experiment, observed for the uncorrelated proton distri- 
bution (dashed line), may be due to the broadening ef- 
fect of temperature on the distribution of the host ions. 
It is worth mentioning that differences between results 
from the present calculations and previous preliminary 
ones 27 are due to refinements of the potentials performed 
presently. 

V. QUANTUM DELOCALIZATION AND THE 
GEOMETRICAL EFFECT 

A. Geometrical Effect vs. Tunneling 

We now address the origin of the huge isotope effect on 
T c observed in KDP and also in the isomorphic H-bonded 
crystals in the family. For forty years, starting from the 
pioneering work of Blinc, 9 the central issue in KDP has 
been whether tunneling is or not at the root of the large 
isotope effect, a fact that was never rigorously confirmed. 
Moreover, a crucial set of experiments pointing against 
the tunneling picture was recently conducted by Nclmcs 
and co-workers: by applying pressure, they conveniently 
tuned the D-shift parameter Sdkdp in DKDP to make 
it coincide with the H-shift parameter Skdp in KDP, 
and they observed that T^ KDP almost coincided with 
Tf DP , in spite of the mass difference between D and H in 
both systems. 18 ' 19 ' 52 This suggests that the modification 
of the H-bond geometry by deuteration - the geometrical 



effect - is a central mechanism in the transition, and is 
intimately connected with the isotope effect. 

FE Phase PE Phase 
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FIG. 6. Momentum distributions along the H-bond in both 
phases of KDP. Experimental data are shown in solid lines. 14 
Left panel (FE phase): calculations for single uncorrelated H 
motion (dashed line), and for H correlated with host relax- 
ation (dotted line). Right panel (PE phase): calculations for 
4-D and 7-D cluster dynamics, as explained in text (dashed 
and dot-dashed lines, respectively). 

As the cluster size grows (N — > oo), the tunnel split- 
ting f2 vanishes. However, only large clusters are ex- 
pected to be relevant for the nearly second-order FE 
transition in these systems. 56 Thus, for large tunneling 
clusters the potential barrier is sufficiently large, and the 
GS levels are deep enough (see Fig. 5) that the rela- 
tion hfljj{D) "C KbT c is satisfied, so much for D as 
for H. Therefore, according to the tunneling model, the 
above relation implies that a simple change of mass upon 
deuteration at fixed potential cannot explain the near 
duplication of T c . In fact, if we consider the largest clus- 
ter (N=10) in Fig. 5 for DKDP, which is larger than 
the crossover length in this system, we have that the GS 
level for the deuterated case (calculated with an effective 
mass of 10/i £> = 35.4 m p ) is around Eqs — -107 meV, 
well below the central barrier, and the tunnel splitting 
amounts to http = 0.34 K. Modifying the mass at fixed 
potential (10/j.h = 25.3 m p ) leads to a tunnel splitting 
of Ml D = 1.74 K. Since Tf KDP « 229if, the relation 
%Q,H(d) *C KbT c clearly holds, and the change in Q, at 
fixed potential accounts only for a small change in T c . 
This is in agreement with the high-pressure experiments 
mentioned above, 18 ' 19 ' 52 where at fixed structural condi- 
tions, the isotope effect in T c appears to be rather mod- 
est. 

Also the geometric effect in the H-bond is very small 



at fixed potential. In fact, the proton and deuteron wave 
functions (WF) in the DKDP potential for the N=7 clus- 
ter, which are reported in Fig. 7 (a), are only slightly dif- 
ferent. As a matter of fact, the distance between peaks 
as a function of the effective mass at fixed potential re- 
mains almost unchanged, as can be seen by the square 
symbols in Fig. 7 (c). 

In contrast, the proton WF for the N=7 cluster in the 
KDP potential exhibits a single broad peak, as shown in 
Fig. 7 (a). But, how can we then explain such a large ge- 
ometric change in going from DKDP to KDP? The first 
observation arises from what is apparent in Figure 5: en- 
ergy barriers in DKDP are much larger than those in 
KDP, implying that quantum effects are significantly re- 
duced in the expanded DKDP lattice. In fact, the proton 
WF in the DKDP potential has more weight in the mid- 
dle of the H-bond (x ps 0) than the deuteron WF (Fig. 
7 (a)). That is, due to zero-point motion, protons are 
more favored in the H-centered position than deuterons. 
This affects the covalency of the bond, which becomes 
stronger as the proton moves to the H-bond center, as 
discussed in Section III. The geometric change of the O- 
H-0 bridge, produced by the combined effect of quan- 
tum derealization and gain in covalency, affects in turn 
the crystal cohesion. Thus, the increased probability of 
the proton to be midway between oxygens strengthens 
the O-H-0 covalent grip and pulls the oxygen atoms to- 
gether, causing a small contraction of the lattice. As will 
be shown in the following Section, this contraction has 
the effect of decreasing the barrier height, thus making 
the proton even more dclocalized. This triggers a further 
contraction of the lattice, and so on in a self-consistent 
way. This self-consistent procedure is finally identified as 
the phenomenon that makes the lattice shrink from the 
larger classical value to the smaller value found for KDP. 
This phenomenon, triggered by tunneling and quantum 
derealization, leads to an enhancement of the geometri- 
cal effect. The overall self-consistent effect is eventually 
much larger than the deuteration effect obtained at fixed 
potential, i.e. at fixed lattice constant. 

To estimate an upper limit to that effect, we compared 
the lattice parameters and the bridge lengths by carrying 
out electronic calculations with classical nuclei (clamped 
nuclei). This was done for two different situations: one 
with the hydrogens forced to stay in the middle of the H- 
bond, and the other with the hydrogens fully off-centered 
in the FE state of KDP. In the latest case, the distance 
between oxygens is doo ~ 2.50 A, falling to doo ~ 2.42 
A when H is centered. In addition, the lattice volume 
is contracted by about 2.3 %. Thus, the proton center- 
ing acts as a very strong attraction center, pulling the 
two oxygens together. We estimate that, at the equilib- 
rium volume, the proton centering creates an equivalent 
pressure of w 20 Kbar. In the true high-temperature PE 
phase, though, the protons are not centered in the mid- 
dle of the H-bonds, but they are equally distributed on 
both sides of the bond, thus reducing the magnitude of 
the effect. 



B. The Isotope Effect: a Nonlinear Self-consistent 
Phenomenon 

In the previous subsection we discussed how a self- 
consistent mechanism combining quantum derealiza- 
tion, the modification of the covalency in the bond, and 
the effect on the lattice parameters, can account for the 
large geometric effect observed upon deuteration. This 
mechanism is now capable of explaining, at least quali- 
tatively, the increase in the order parameter and T c with 
deuteration. This self-consistent mechanism has obvi- 
ously its origin in the difference in tunneling induced by 
different masses, but is largely amplified through the ge- 
ometric modification of bond lengths and energy scales. 

To demonstrate the effect of isotopic substitution 
via this self-consistent non- linear mechanism, we con- 
structed the following simple model: we considered 
the Schrodingcr equation for the clusters with a WF- 
dependent term added to the bare potential. The effec- 
tive potential reads: 



V eS (x) = V (x)-k\*(x)\' 



(2) 



where x is the collective coordinate of the cluster and 
Vq(x) is a quartic double- well similar to those of Fig. 5. 
The term in ^(rc)! 2 serves as a non-linear feedback in 
the model: when the particle is more delocalized, it has 
more weight in the middle of the H-bond. Then, ^(x)! 2 
increases at the center, the effective barrier is lowered, the 
particle further delocalizcs, and so on, self-consistently. 
The bare potential can be written as 



V (x) 



2x 



+ 



2x 



(3) 



in terms of its energy barrier E® and minima separation 
S^in- The parameters values k = 20.2 meV.A, E® = 
35 meV and 5° min = 0.24 A were chosen so as to quali- 
tatively reproduce the WF profiles in the cases of KDP 
(broad single peak) and DKDP (double peak), for the 
same cluster size. Once these parameters are fixed, the 
WF self-consistent solutions depend only on the effec- 
tive mass. Figure 7(b) shows the WF corresponding to 
yU_D (solid line) and Hh (dashed line), which are similar 
to those calculated from the ab initio potentials for the 
N=7 cluster (Fig. 7(a)). 

In Fig. 7(c), we show the distance between peaks 5 P 
in the WF as a function of the cluster effective mass /z. 
Starting from the finite value for no (DKDP), S p de- 
creases remarkably towards lower [i values, until it van- 
ishes near /Ujj (KDP) (see circles in Fig. 7(c)). This 
strong dependence of 5 P on the mass is in striking con- 
trast with the very weak dependence obtained at fixed 
DKDP potential and geometry (square symbols). Such a 
large mass dependence, can now explain the large isotope 
effect found in KDP, via an amplified and self-consistent 
geometrical modification of the H-bond. 
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FIG. 7. WF in the 7-H(D) cluster PES for (a) a& initio and 
(b) self-consistent model calculations. Solid (dashed) lines are 
for D (H). Dotted line is for H in the DKDP PES. (c) WF 
peak separation S p as a function of the cluster effective mass 
fj, (given in units of the proton mass) for the self-consistent 
model (circles) and for fixed DKDP potential (squares). Lines 
are guides to the eye. 



VI. PRESSURE EFFECTS 



culations and the aid of the previously introduced self- 
consistent model for the geometrical effect. 53 To this pur- 
pose, we first focus on first-principles calculations, where 
we used the PW approach explained in section II, com- 
bining ultra-soft (O and H) and norm-conserving (K and 
P) pseudopotentials. 54 ' 39 Using linear-response theory, 40 
the unstable ferroelectric mode at the zone-centre was 
identified. This corresponds mostly to H-atoms displace- 
ments with the pattern shown in Fig. 3, with heavy ions 
displacing to a lesser extent. The O-atoms are practically 
fixed in this mode. 

The FE mode amplitude is identified with the H off- 
centering coordinate (x). The total energy profile as a 
function of x displays an effective double-well potential 
for the FE mode. 46 Thus, the potential can be character- 
ized by two parameters: the energy barrier Ef, between 
the stable and the unstable (x = 0) configurations, and 
the separation between minima S m . The values of Ef, vs. 
5 m obtained under different applied pressures are plotted 
in Fig. 8. A nearly quadratic behaviour with simultane- 
ous vanishing of Ef, and 5 m is observed. Classically, ferro- 
electricity would dissapear above w 100 Kbar. However, 
the critical temperatures vanish at substantially lower 
critical pressures P c , which are isotope- dependent: 17 
Kbar for KDP and 60 Kbar for DKDP. 55 This can be 
understood by considering the quantum character of the 
nuclear dynamics. In fact, the zero-point energy, which 
is larger for the proton, should lower the effective energy 
barrier leading to lower critical pressures. 



Experiments under pressure carried out during the late 
eighties showed, within error bars, that T c depends lin- 
early on 5 for different H-bonded ferroelectric materi- 
als. 52 Strictly speaking, the linear relation is verified for 
some of these compounds within a restricted region of the 
T c vs. S plot. In the case of KDP, the linear behaviour 
appears to extend up to a pressure of 17 Kbar, where 
T c vanishes. The need for very high pressures (w 60 
Kbar) to achieve a vanishing T c for DKDP and other re- 
lated H-bonded compounds, precluded the generalization 
of the enunciated hypothesis, and motivated the devel- 
opment of improved diamond anvil cells. Nevertheless, 
a striking observation arises when the extrapolation of 
the linear behaviour in the mentioned materials is car- 
ried down towards lower values of T c : the critical tem- 
perature appears to vanish for all systems, deuterated 
and protonated, with one, two and three-dimensional H- 
bond networks, around a seemingly universal point where 
5 = S c « 0.2 A. 52 

In particular, the effect of deuteration can be reverted 
by applying pressure, and the critical temperature of 
KDP can be reproduced by compressing DKDP in such 
a way that the structural parameter Sdkdp assumes a 
value very close to that measured in KDP at the initial 
pressure. 52 This is valid at all the measured pressures. 

In this Section we explore the connection between pres- 
sure and isotope effect by means of first-principles cal- 
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FIG. 8. Energy Barrier E& as a function of the distance 
between the two minima of the double- well, S m , for different 
pressures. The dashed line is a guide to the eye. An approx- 
imate conversion to temperatures has been included to the 
right of the graph. 



In Fig. 9 we show the two parameters of the double- 
well, i.e. Ef, and S m , as a function of pressure. In this fig- 
ure, critical pressures for KDP and DKDP are indicated 
by vertical arrows. These correspond, in our calculation, 
to different classical values of the minima separation S m , 
which are around 0.3 A in KDP and 0.18 A in DKDP. 
Experiments, however, indicate that 6 should approach 
the same universal value 5 C ~ 0.2 A for both compounds 
as T c goes to zero. The difference found in the S m values 
should again compensate for the quantum correction due 
to the nuclear dynamics. 




100 120 



Pressure [Kbar] 



FIG. 9. Energy barrier Ef, (triangles) and double-well min- 
ima separation 5 m (squares) as a function of pressure. Lines 
are guide to the eye only. 

In fact, being lighter, protons will localize closer to the 
middle of the H-bond than deuteriums, leading in princi- 
ple to such compensation. However, we will show below 
that this effect alone is not sufficient. The self-consistent 
geometrical effect discussed in the previous section, which 
was essential to explain the huge variation in the order 
parameter upon deuteration 29 is also crucial to explain 
the close similarity of 8 C for deuterated and protonated 
systems. 

Neutron diffraction experiments indicate that, to have 
the same value of 6, different pressures have to be ap- 
plied to KDP and DKDP. 52 When converting pressures 
into global energy barriers using Fig. 9, the difference in 
energy barrier required to have the same value of 5 turns 
out to be nearly independent of S, assuming a value of 
ss 23 meV per unit cell (two formula units) . This energy 
difference then seems to be a key quantity, which takes 
into account that the WF of the more-easily-tunneling 
protons in KDP will exhibit the same distance between 
peaks as the WF for deuterons in DKDP, only if the un- 
derlying double-well is significantly deeper, i.e. by 23 



meV, and also the distance between minima S m is in- 
creased. The reason for this increased separation is in 
the very nature of the hydrogen bond: the more distant 
the O-atoms, the less covalent the bond, the larger S m , 
and the deeper the double- well. This is the essence of the 
geometrical effect, and 23 meV is the energy difference 
required to modify the 0-0 distance in such a way that 
the distance between peaks in the WF is the same for 
KDP and DKDP. 

To show how the geometrical effect enters into play, we 
considered the feedback effective potential V e s(x) from 
Eq. 2. The probability distribution for the H(D) motion 
^(x)! 2 , was obtained by solving the Schrodinger equa- 
tion in the effective potential V e g(x). We chose the N=7 
cluster, with the relation hd/hh ~ 1-3 for the effective 
masses in DKDP and KDP. The value of 5H, was fixed to 
. Within this model, we studied how the value of 
the energy barrier E® has to be modified (simulating the 
application of pressure), in order to keep the peak sepa- 
ration 8 of the wave function ^>(x) constant upon deuter- 
ation. This study was carried out for different values of 
the non- linear parameter k in the model. Large values 
of k represent important feedback geometric effects. 
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FIG. 10. Average pressure P as a function of the coupling 
constant k from the non-linear model, as defined in text. Lines 
are guide to the eye only. 

For each value of the non-linear parameter k, we 
searched for values of the barrier energies E^ and E® 
for each isotope such that S KDP = S DKDP , with the ad- 
ditional constraint that the difference E^ — E® = 23 
meV remains constant. Using Fig. 9, we converted en- 
ergy barriers into pressure, with the warning that these 
correspond to global energies, while E° represents clus- 
ter energy barriers. Since here we are interested in 
qualitative issues, the former is a reasonable approxi- 
mation. Therefore, we calculate the average pressure 
P = {P(Ef) + P{Ef?)} /2 for each value of k, which 



serves to fix the absolute scale of pressures necessary to 
fulfill the above constraints. In Fig. 10 we plot P as a 
function of k. For a value of k = (no geometrical ef- 
fect), in order to maintain S constant upon deuteration, 
we have to apply a large negative pressure to the system, 
and hence expand it to have significantly higher energy 
barriers. Conversely, for larger pressures, as those exper- 
imentally measured, compensation can be achieved only 
by considering large values of k, thus leading to impor- 
tant structural non-linear effects. In summary, to explain 
pressure effects on KDP and DKDP, 52 it is necessary to 
consider the non-linear relation between isotope substi- 
tution and geometric effect. 

VII. DISCUSSION AND CONCLUSIONS 

The nature of the instability that drives the FE tran- 
sition and leads to the onset of the spontaneous polariza- 
tion P s has been extensively discussed in the past. 10 ' 57 ' 45 
Although the H(D) ordering in the basal plane is un- 
doubtedly correlated with the transition, it was origi- 
nally assumed that P s , which is oriented along the c-axis, 
was due to the displacements of K + and P~ ions along 
this axis. 10 However, the observed value of P s can only 
be explained within this model if unrealistic, very large 
charges for the phosphorous ion are assumed. 45 Bystrov 
and Popova 45 proposed, alternatively, that the source of 
P s could be the electron density shift in the P-0 and P- 
O-H bonds in the polar direction, which occurs when the 
protons order perpendicularly. This assumption cannot 
be assessed through model calculations, for it is origi- 
nated in the complex electronic interactions in the sys- 
tem. 

By means of the present ab initio calculations we were 
able to overcome this limitation, and to show that the 
FE instability has its origin on an electronic charge re- 
organization within the internal P-0 and P-O-H bonds 
of the phosphates, as the H-atoms order off-center in 
the H-bonds. As a matter of fact, the overall effect 
produced by the H-ordering is an electronic charge flow 
from the 02 side to the 01 side of the PO4 tetrahedron, 
and a concomitant distortion of the former. 26 This is in 
agreement with the explanation given in Ref. 45 , and also 
agrees with the results obtained by another recent first- 
principles calculation. 28 

The microscopic origin of the global FE instability, i.e. 
the connection between H ordering and phosphate dis- 
tortions, is also demonstrated in the strong correlation 
between the off-centering parameters S for H and 7 for 
K-P (see inset of Fig. 4). One of the observations here 
is that the overall potential for the H motion is not in 
fact separable, and one has to deal with the problem of 
the "chicken and the egg", what is really first? 22 Nev- 
ertheless, using that 7 is a measure of polarization, we 
showed that the source of the FE instability is in the H 
off-centering, and not vice versa. 



There is a long controversy about the origin of the 
FE transition. Some experimental facts support the 
coupled proton-phonon model which displays essentially 
a displacive-like transition. 58 Other experiments, e.g. 
Raman studies 59 , seem to indicate the importance of 
the order-disorder character of the transition originated 
in the H2PO4 unit dipoles. Electron- nuclear double- 
resonance (ENDOR) measurements 60 indicate that not 
only the H2PO4 group, but also the K atoms, are disor- 
dered over at least two configurations in the paraelectric 
phase. It is also shown by neutron scattering experi- 
ments that the P atom is distributed over at least two 
sites in DKDP. 50 ' 51 In spite of the still unresolved char- 
acter of the transition, it is clear that local instabilities 
arising from the coupling of light and heavy ions are very 
important in this system, irrespective of the correlation 
length scale associated with the transition. 

In fact, our calculations in the PE phase show that 
local proton distortions with the FE mode pattern need 
to be accompanied by heavy ion relaxations in the PO4- 
K group to produce significant instabilities, a fact which 
is in agreement with experiments. We have shown that 
the correlation length associated with the FE instability 
is much larger in KDP than in DKDP, suggesting that 
DKDP will behave more as an order-disorder ferroelectric 
than KDP. 

The coherent interference of the proton in the two 
equivalent sites of the PE phase was observed recently 
by neutron Compton scattering experiments. 14 Quan- 
tum coherence arises in our calculations for DKDP only 
when P and K ions are allowed to relax together with the 
deuterons. In KDP, the onset of tunneling and hence, co- 
herence, would require relaxations of clusters comprising 
more than three KH2PO4 groups, which were not consid- 
ered in the present work. The momentum distributions 
calculated in the PE phase are in qualitative agreement 
with the experiment. 14 Quantum coherence would, thus, 
be produced by a dressed proton, i.e. strongly correlated 
with the heavier ions. In fact, in Ref. 14 , many-body ef- 
fects due to the motion of the surrounding ions are not 
excluded, but it turns out that these are difficult to as- 
sess. We have also found good agreement with the exper- 
iment for the momentum distribution in the FE phase, 
which corresponds to a single, anharmonic well for each 
individual proton, in a host FE lattice. In this situation 
there is no coherence between the motions of the various 
protons. 

The most striking feature, which is not yet satisfac- 
torily understood, is undoubtedly the huge isotope ef- 
fect in the critical temperature and the order parameter 
of the transition. The first explanation was that pro- 
posed by the tunneling model and later modifications, 9 ' 10 
but soon after the vast set of experiments carried out 
by Nelmes and co-workers, 47 > 18 > 19 > 50 > 52 and the compre- 
hensive structural compilation of Ichikawa et at, 16 the 
importance of the so- called geometrical effect as an al- 
ternative explanation became apparent. Other experi- 
ments 59,15 ' 14 and models 11 ' 23 ' 13 favoured one or the other 



vision, or even both, but an overall and consistent expla- 
nation of the phenomenon is still lacking. Still unan- 
swered questions like: if tunneling occurs, what are the 
main units that tunnel?, what is the connection between 
tunneling and geometrical effects?, and what is the true 
microscopic origin of the latter?, are possibly some of the 
reasons why a full explanation of the isotope effect is not 
yet available. With the aid of the ab initio scheme, our 
efforts in this direction shed also light onto the underly- 
ing microscopic mechanism for the isotope effect. 

Protons alone are not able to tunnel because the ef- 
fective potentials in which they move display tiny double 
wells, and the particles are broadly delocalized around 
the center of the H-bonds. In this regard, we conclude 
that the simplified version of the tunneling model, i.e. 
that of a tunneling proton, or even a collective proton 
soft-mode alone, is not supported by our calculations. On 
the other hand, we observe "tunneling clusters", with an 
effective mass much larger than that of a single, or even 
several, protons (deuterons), due to the correlation with 
heavier ions. These clusters have different sizes, leading 
to different lengths and energy scales competing in the 
system in the PE phase. The smallest tunneling unit 
in DKDP is found to be the KD2PO4 group. This re- 
sults agree with the idea developed by Blinc and Zeks 
of a tunneling model for the whole H2PO4 unit, which 
helped to describe the typical order-disorder phenomena 
observed in some experimental trends. 58 However, the 
explanation for the isotope effect arising in the present 
work is even more complex, and goes beyond the concept 
of "tunneling alone" as its main cause. 

Although the PE phase of the system shows a complex 
scenario due to the appearance of different length scales, 
it is clear that larger clusters will prevail as the transition 
approaches. We have shown for both isotopes that tun- 
nel splittings in these clusters at fixed potential are much 
smaller than the thermal energy at the critical tempera- 
ture. Thus, at fixed potential, tunneling is not able to ac- 
count for the large isotope effect in the system. However, 
as the dressed particle is delocalized by tunneling, the 
effective potential felt by the H(D)-atom changes upon 
isotopic substitution, due to significant modifications in 
the chemical properties of the 0-H...0 bond, which are 
reflected in a concomitant lattice relaxation. 

With the aid of a simple model based in our ab initio 
results, we were able to show how this feedback effect 
strongly amplifies the geometrical modifications in the 
H(D)-bridge. Tunneling triggers a self-consistent mecha- 
nism, but in the end, the geometrical effect dominates 
the scenario and accounts for the huge isotope effect, 
in agreement with neutron scattering experiments. 19 ' 52 
Therefore, these aspects, which were largely debated in 
the past, here appear as complementary and deeply con- 
nected to each other. 11 ' 57 ' 29 

The feedback effect of the geometrical modifications on 
the proton distribution is also necessary to explain the 
results of experiments under pressure. 52 There, it was 
observed that the critical pressures at which the transi- 



tion temperature vanishes correspond to an isotope and 
material-independent value of the peak separation in the 
proton distribution. We have shown that this unique 
value for KDP and DKDP can be achieved only as a 
consequence of a compensation between quantum der- 
ealization effects and geometrical modifications imposed 
by pressure. 53 

The question of why T c is so closely related to the 
distance between H peaks 8 is still open. A possible ex- 
planation can rely on the fact that equal GS levels rela- 
tive to the top of the barrier in a double well potential 
correspond to approximately equal 8, irrespective of the 
energy barrier height and mass of the tunneling particle. 
We have verified in test calculations that this holds as 
long as the GS energy is not very deep below the top of 
the energy barrier, which means not too low pressures, as 
discussed in Section VI. In addition, in our mean-field de- 
scription of the scenario of the FE transition, T c should 
be related to the energy difference between the GS level 
and the top of the barrier in the double well potential 
of a tunneling cluster (see Fig. 5). This could explain 
the close relation between T c and 5 observed in neutron 
diffraction experiments. 52 

The nonlinear feedback between tunneling and struc- 
tural modifications is a phenomenon of wider implica- 
tions. Tunneling units are indeed observed in a large 
variety of molecular compounds and biomolecules. Both 
tunneling and structural changes are important for the 
reaction mechanisms of enzymes 62 and other biological 
processes. Our results on KDP supports the already ex- 
pressed need for revision of the general theories of host- 
and-tunneling systems. 22 

In summary, we showed that proton ordering in KDP 
leads to an electronic charge redistribution and ionic dis- 
placements that originate the spontaneous polarization 
of the ferroelectric phase. The instability process is con- 
trolled by the hydrogen off-centering. The double-peaked 
proton distribution in the bridges, observed in the para- 
electric phase, cannot be explained by a dynamics of pro- 
tons alone. These must be correlated with displacements 
of the heavier ions within clusters. These tunneling clus- 
ters can explain the recent evidence of tunneling obtained 
from Compton scattering measurements. We also showed 
that the mere mass change upon deuteration docs not 
explain the huge isotope effect observed. We find that 
structural changes arising from the modification of the 
covalency in the bridges produce a feedback effect on the 
tunneling that strongly enhances the phenomenon. The 
resulting influence of the geometric changes on the iso- 
tope effect is in agreement with experimental data from 
neutron scattering. Moreover, the behavior of the pro- 
ton/dcutcron distribution in the bridges under pressure 
can only be explained by invoking the mentioned feed- 
back effect of geometry. 
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